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Summary. The current use of functionals to evaluate order-of-convergence of a nu- 
merical scheme can lead to incorrect values. The problem comes about because of 
interplay between the errors from the evaluation of the functional, e.g., quadrature 
error, and from the numerical scheme discretization. Alternative procedures for de- 
ducing the order-property of a scheme are presented. The problem is studied within 
the context of the inviscid supersonic flow over a blunt body; however, the problem 
and solutions presented are not unique to this example. 


1 Introduction 

Computational Aerodynamicists conduct most of their grid convergence stud- 
ies by studying the behavior of solution functionals, e.g., drag, lift and moment 
coefficients, as the computational grids are refined. Functionals are used for 
several reasons: first, their accurate evaluation is of intrinsic value; and second, 
they provide a means of determining convergence properties of a numerical 
scheme without looking directly at hundreds of thousands of field point val- 
ues. Ideally, an error measure should be used to examine order-properties of 
grid convergence studies; however, exact solutions are usually not available for 
flows of practical interest. Therefore, estimating convergence properties using 
functionals is frequently the only course of action available. 

However, there are some subtle problems associated with the use of func- 
tionals for grid convergence studies, and if these problems are not recognized 
and resolved, the results that follow from the use of functionals can be very 
misleading. It is the purpose of this paper to expose these problems, and where 
possible, suggest solutions. 

There are many aspects of a numerical order-properties analysis that must 
be done correctly in order for the analysis to be reliable. Paramount among 
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these are: that grid refinements must be uniform, preferably with grids se- 
quences that are hierarchical; and the iterative methods used to solve the 
discrete equations on a given grid must be sufficiently converged, preferably 
with residuals reduced several orders of magnitude below the solution error. 
Of course, this is complicated by the fact that the errors are not known a 
priori. 

The problems that are associated with the use of functionals for the study 
of grid convergence rates are illustrated with numerical results from the com- 
putation of a blunt-body in an inviscid supersonic stream. However, it should 
be emphasized that the problems discussed are not unique to the blunt-body 
problem, indeed they are not unique to fluid mechanics, and may occur in any 
grid convergence study involving functionals. The particular case studied is 
the Mach 6 flow of an inviscid gas over a circular cylinder. In the numerical 
implementation the problem is solved as a time dependent problem with the 
bow shock wave fitted as a boundary of the flow. By fitting the shock, the 
numerical scheme acts only on a smooth flow region. Thus, the computation 
is limited to the layer bounded by the bow shock, the circular cylinder, the 
symmetry line, and a supersonic outflow boundary imposed at some 9 = 0 max , 
see Fig. 1. 



x 


Fig. 1 . Supersonic blunt-body flow field, showing isobars, Moo = 6. 


The physical plane is transformed to a computational plane where N and 
M mesh points are uniformly distributed between the body and the shock 
and between the symmetry line and the outflow line, respectively. The predic- 
tor/corrector MacCormack scheme [1] is used for the numerical integration of 
the Euler equations. 

Table 1 shows results obtained with a series of grids. The number of mesh 
points corresponding to each k th grid are N*, = 3 x 2 fc and M& = 5 x 2 k . 
Columns 2 and 3 display the inviscid drag coefficient computed with the 
trapezoidal rule (TR) and with Simpson’s rule (SR). The drag coefficient 




Problems Associated with Grid Convergence of Functionals 


3 


orcler-of-convergence is given by 


Pk = log 2 


C d,k C d^k + 1 

Cd ,fc+ 1 - Cd ,/c+2 J 


(i) 


and is shown in the last two columns. The order-of convergence for k = 3 
for TR and SR shows a large discrepancy and both results are significantly 
greater than the formal order of the scheme which is second order. For k = 4, 
the drag coefficient is not monotone and the order-of-convergence evaluation 
fails. (Note that the order-of-convergence for grid k depends on the solutions 
from grids k, k + 1, and k + 2). 


Table 1. M ^ = 6 cases investigated. For the three finest grids the trapezoidal rule 
(TR) computed drag is not monotone and for both quadrature rules the computed 
drag exhibits super-convergence. 


k 

C d (TR) 

C d (SR) 

P(TR) p(SR) 

1 

1.8755919 

1.8767669 

1.92 

1.90 

2 

1.8706109 

1.8709412 

2.73 

2.58 

3 

1.8692925 

1.8693766 

4.03 

3.25 

4 

1.8690942 

1.8691147 

— 

— 

5 

1.8690821 

1.8690872 



6 

1.8690859 

1.8690872 




To establish that there is reason to suspect these results, consider the 
behavior of the error norm in total temperature. For this problem in the steady 
state, the total enthalpy, and hence the total temperature, is constant. The L 2 
and Lqo norms of the total temperature error are shown in Fig. 2. The order- 
of-convergence based on the L 2 and L^ norms is 2.03 and 1.84, respectively. 
These are in fairly good agreement with the formal order of the scheme. Why 
then is the order-of-convergence of the drag functional misbehaving? 


2 Order-of-convergence of functionals 

To answer the last question, first consider the following question: If the sur- 
face pressure converges with order p , what should the expected order-of- 
convergence of the drag-functional be? To this end, let the computed sur- 
face pressure, P c , normalized by the free stream pressure, Poo, be given by 
Pc/Poo = P&! Poo + a(6)h p , where P e is the exact surface pressure. The sec- 
tional drag coefficient is defined by 

f^max 

c d = / (p c /p 00 -i)b(e)cos(e)de/(^MlA ref ) ( 2 ) 

Jo 
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Fig. 2. Z/2 and L x of total temperature error for entire shock layer, based on results 
from grids k = 3,4,5, with hk = 1/y/NkMk- 

where 7 is the ratio of specific heats, is the free stream Mach number, b(9) 
is the body radius, and A re f is a reference surface area, here taken as the pro- 
jected plan- form area. Therefore, Cd = Cd, e + rh p where r is a constant and 
Cd, e is the exact value of the drag coefficient; therefore, the drag-functional 
should converge with order p. 


3 The problem with quadrature 

If it is assumed that the pressure order-of-convergence behaves like the total 
temperature order-of-convergence, then the result just obtained for drag is 
not consistent with the results of Table 1. The problem lies in the numerical 
integration of (2) . The integration is approximated by a quadrature taken over 
M equally spaced intervals on the surface of the cylinder, i.e. f Q max fdO ss 
aifi . The quadrature has a leading error of order h q . It is easy to show 
that Cd = Cd, e +Pqh q + rh p + 0{h q+p ), where [3 q is a constant and q equals 
2 for TR and 4 for SR. Using this relation for Cd in (1) we find 

__ f 2 p ~ 2 /3 q [1 — 2 q ] + h p ~ q T [1 — 2 P ] \ 

P P + og 2 | 22(p _ 2) ^ _ 29] + h p - q r [1 _ 2P] j ’ ^ 

where h is the coarse grid spacing. Here p is the computed order-of-convergence 
and the log 2 -term is an error brought about by the interplay between the 
quadrature error and algorithmic error. In the limit h — > 0, the p behavior is 
given by p — » q if p > q, and p — > p if p < q. However, in a computation h will 
always be finite and having p < q is not a guarantee that the log 2 -term will 
be small. 
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4 How to eliminate the quadrature error 

A solution to the quadrature problem can be found by studying (1). Consider 
the numerator. The numerator is the difference between the drag coefficients 
of the medium (k+1) and coarse (k) grids. The medium grid has a quadrature 
error of order (h/ 2) 2 , while the coarse grid has a quadrature error of order h 2 . 
These errors do not cancel out and their interplay with the algorithmic error 
causes some (not all) of the problems in the results of Table 1. The solution is 
to implement the quadrature in such a way that the quadrature errors of the 
medium and coarse grids cancel. To do this, evaluate both the medium and 
coarse grid quadratures using an h interval, i.e., use only every other point of 
the medium grid. The same idea is applied to the denominator by evaluating 
both the medium and fine grid quadratures using an h/2 interval. 


5 Higher-order algorithmic error model 


The actual algorithmic error in any numerical solution on any given grid 
contains a full hierarchy of errors that are unknown, but are generally as- 
sumed to be of the form u Ct k = u e + ^r5°=p a «^fc’ where p is the unknown 
actual order of the numerical method. The standard method for deducing or- 
der properties from grid convergence, described earlier, is obtained by fitting 
a single error mode of the form u Ct k = u e + ah p to the actual error. For suf- 
ficiently small h, the actual error is dominated by the lowest order term, and 
the single mode model provides a good fit and an accurate prediction of the 
order-of-convergence. However, for larger h above this asymptotic region (it is 
surprising how small h has to be to reach asymptotic convergence), multiple 
error modes are competing, and their projection onto a single mode can be 
erroneous. Consider then a two- mode error model: 


Cd,k=C d:ex + a 1 h p k + a 2 h p k +1 . (4) 

For this model, using a sequence of four grids (k through k-3), where 
hk/hk-i = hk-i/hk -2 = hk- 2 /hk -3 = 2, we find the order-of-convergence 
to be 


P = log 2 


32\ fc _ 2 , fc-i — w 9 Z \|_ 2 k _ 1 + 


k-2,k-l 

4Z\fc-i,fc 


J fc-l,fe^fe-2,fc-3 


(5) 


where A it j = — C ( ij . It is important to monitor the ratio a 2 h/ai. 

a 2 h _ 4(1-2 p) (2M,_ 2 , t _ 1 -4- 1 ,t) ,, 

Q'l (1 — 2 P+ 1 ) (Ak-2,k — t 2‘ P+l Ak-2,k-l) 

When this ratio is less than one, the one-mode error model is valid. With 
a two-mode error model and eliminating the quadrature error as previously 
described, we obtain the results listed in Table 2. For more details see [2]. 
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Table 2. Drag order-of-convergence using higher order method and coarse grid 
interval, k = 3, for quadrature rules. 


k 

C d (TR) 

C d (SR) 

P(TR) p(SR) 

3 

1.8692925 

1.8693766 

1.674 1.681 

4 

1.8690328 

1.8691145 


5” 

1.8690056 

1.8690872 


IT 

1.8690056 

1.8690871 



6 Conclusions and recommendations 

With the increased reliance in both science and engineering on the numerical 
solution of partial differential equations, the subject of code verification has 
become increasingly significant. An important element of code verification is 
the study of grid convergence. Most studies today of this subject have been 
at best superficial and in many cases painfully inadequate. This paper is an 
attempt to reverse this trend by first highlighting a series of problems that ex- 
ist in the standard order-of-convergence analysis, particularly as it relates to 
the evaluation of functionals, and second by providing a number of solutions 
and workarounds to these problems. It is important to distinguish between 
a code verification effort and an effort to determine if a particular solution 
to a specific problem is sufficiently accurate for some intended use. The two 
tasks are very different. A rigorous grid convergence and order-of-convergence 
study can aid in determining if an algorithm has been implemented correctly. 
However, such a rigorous study requires grids of the same family and grid 
refinements that are uniform, preferably with grids sequences that are nested. 
In the second task, limited time and resources often lead to compromising one 
or more attributes of a rigorous study. While non-uniform mesh refinement 
may lead to some improvement in the solution, order-of-convergence prop- 
erties computed from non-uniform refinements or ill-converged solutions sets 
are meaningless. Whenever possible, error norms should be used to establish 
the order-of-convergence. The higher order analysis developed in Sect. 5 is the 
best way to evaluate if the asymptotic range has been reached or if more levels 
of grid refinement are needed to reach it. It should be part of any rigorous 
grid convergence study. Reference [2] provides a more in depth study of these 
issues. 
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